home *** CD-ROM | disk | FTP | other *** search
/ PC-Blue - MS DOS Public Domain Library / PC-Blue MS-DOS Public Domain Library - NYACC.iso / vol270 / laplace.c < prev    next >
Encoding:
C/C++ Source or Header  |  1986-12-16  |  26.3 KB  |  953 lines
  1. #define SWITCH
  2. /*    laplace - solve Laplace's equation for specified boundary conditions
  3.                 in two dimensions
  4.  
  5.     history...
  6.         28 Jun 86    -g switch sets edge of grid to zero potential, 
  7.                     -m switch specifies a margin around the specified boundaries.
  8.         21 Jun 86    Allowing boundary to be between grid points.
  9.         31 May 86    showing boundries and calculated potentials are optional.
  10.  
  11.     approximate performance...
  12.  
  13. h=w=32, n=64
  14. Current time is 20:58:25.76
  15.  
  16. Current time is 21:00:03.36
  17.  
  18.  
  19. h=w=64, n=128
  20. Current time is 21:01:11.29
  21.  
  22. Current time is 21:07:08.87
  23.  
  24. ...producing 42K file
  25.  
  26.         ...on Z-100 with V20 at 7.5 MHz, all files on ramdisk.
  27.  
  28.     notes...
  29.         need to update above performance.
  30.         should investigate recursive subdivision by factors other than 2.
  31.         needs -m switch to specify a margin (for expanding the x & y range)
  32.         can't handle boundary crossing edge of grid: referencing nonexistant
  33.             grid points.
  34.         needs option for cylindrical symmetry
  35.  
  36. */
  37. #define DESMET            /* deleting this only disables the timing function */
  38.  
  39. #include <stdio.h>
  40. #include <math.h>
  41.  
  42. #define VERSION "0.2"
  43. #define MAX_CROSS 300        /* max # times boundary can cross grid */
  44. #define ENTRIES 300            /* maximum # points in boundary point curves */
  45. #define MAXBOUNDARIES 50    /* maximum # separate curves in boundary */
  46. #define PERMITTIVITY 8.85418e-12    /* permittivity of free space, farad/m */
  47. #define SCALE 256
  48. #define MAX(a,b) ((a)>(b)?(a):(b))
  49.  
  50. typedef struct
  51.     {int grid;    /*    cell number */
  52.     int direct;        /*    direction: 0, 1, 2, 3  for  +x, -x, +y, -y  */
  53.     double dist;    /*    distance from cell center and boundary crossing point 
  54.                         (fraction of grid spacing)  */
  55.     int volt;        /*    negative of boundary # (so it's an index into volt[]) */
  56.     int dummy[4];    /*    ...so sizeof(crossing) >= sizeof(boundary_rule)    */
  57.     } crossing;
  58.  
  59. typedef struct
  60.     {int var[4];    /*    grid number on which this point depends    */
  61.     int coef[4];    /*    corresponding coefficient, scaled up by SCALE */
  62.     } boundary_rule;
  63.  
  64. union     /*    these can share space */
  65.     {boundary_rule    b_rule[MAX_CROSS];
  66.     crossing        cross_list[MAX_CROSS];
  67.     } bo;
  68.  
  69. show_b(b) boundary_rule *b;
  70. {    int i, sum;
  71.     sum=0;
  72.     for (i=0; i<4; i++) sum += b->coef[i];
  73.     printf("rule %d: <%2d>*%d+<%2d>*%d+<%2d>*%d+<%2d>*%d  (ttl wt=%d)\n",
  74.         b-bo.b_rule,
  75.         b->var[0],b->coef[0],
  76.         b->var[1],b->coef[1],
  77.         b->var[2],b->coef[2],
  78.         b->var[3],b->coef[3],
  79.         sum);
  80.     i=0;
  81. }
  82.  
  83. char *c_label[4]={"right of","left of","above","below"};
  84. show_c(c) crossing *c;
  85. {    printf("boundary value %d is %f %s grid %d\n",
  86.     c->volt, c->dist, c_label[c->direct], c->grid);
  87. }
  88.  
  89.  
  90. int lc;            /*    index of 1st unused entry in bo.cross_list    */
  91. int lb;            /*    index of 1st unused entry in bo.b_rule    */
  92.  
  93. /*
  94.     cell types are as follows...
  95.  
  96.             2 3 3 3 4
  97.             9 1 1 1 5
  98.             9 1 1 1 5
  99.             9 1 1 1 5
  100.             8 7 7 7 6
  101.  
  102.     ...and type 0 cells have a specified potential (part of the
  103.     boundary conditions) and are never updated.  There are special cell
  104.     types are for boundaries where the normal gradient vanishes:
  105.  
  106.              10  10  10  10
  107.            13              11
  108.            13              11
  109.            13              11
  110.              12  12  12  12
  111. */
  112. int alpha=0, alphap1=1;
  113. int *type; int *volt;
  114. int debugging=0;
  115. int show_resid=0;    /* nonzero if residuals are to be shown */
  116. int grounded=0;        /* nonzero if edge of grid is to be grounded */
  117. int show_boundaries=0;    /* nonzero if boundaries are to be shown */
  118. int show_potentials=0;    /* nonzero if calculated potentials are to be shown */
  119. int left_symmetric=0;    /* nonzero if x gradient vanishes on left edge of grid */
  120. int lower_symmetric=0;    /* nonzero if y gradient vanishes on lower edge of grid */
  121. int done=0;
  122. long work=0;        /* number of cell updates so far */
  123.  
  124. double xmin=1.e30;
  125. double xmax=-1.e30;
  126. double ymin=1.e30;
  127. double ymax=-1.e30;
  128. double zmin=1.e30;
  129. double zmax=-1.e30;
  130. double margin=0.;        /* minimum amount to expand grid beyond given boundaries */
  131. double *x;                /* pointer to data area for boundaries */
  132. double *y;                /* pointer to data area for boundaries */
  133. double xscale,yscale,zscale;
  134. double edge;
  135. double aspect,correct;
  136. double p_voltage[MAXBOUNDARIES];
  137. double residual();
  138.  
  139. int boundaries=0;        /* number of boundaries in input data */
  140. int x_arguments=0;
  141. int y_arguments=0;
  142. int height=32;        /* default height of grid in cells */
  143. int width=32;        /* default width of grid in cells */
  144. int cycles=64;        /* default # relaxation cycles */
  145. int n;                /* number of entries in x and y */
  146. int index_array[MAXBOUNDARIES];    /* indexes into x and y */
  147. int *p_data=index_array;
  148. char outname[40];
  149.  
  150. FILE file;
  151. FILE ofile;
  152.  
  153.  
  154. main (argc,argv) int argc; char **argv;
  155. {    int i,j,k,hw,col;
  156.     double yval;
  157.     char *s;
  158.     read_data(argc,argv);
  159. /*
  160.     printf("there are %d points...\n",n);
  161.     k=0;
  162.     for (i=0; i<n; i++)
  163.         {printf("%f %f\n",x[i],y[i]);
  164.         if(i==p_data[k]) printf("       ...at voltage %f\n",p_voltage[k++]);
  165.         }
  166. */
  167.     if(s=index(argv[1],'.')) strncpy(outname,argv[1],s-argv[1]);
  168.     else strcpy(outname,argv[1]);
  169.     strcat(outname,".pot");
  170.     ofile=fopen(outname,"w");
  171.     if(ofile==0) {printf("can\'t create output file\n"); exit(1);}
  172.     aspect=(ymax-ymin)/(xmax-xmin);
  173.     if(height<width*aspect) height=width*aspect;
  174.     else width=height/aspect;
  175.     if(height<width*aspect)
  176.         {correct=(ymax-ymin)*width/(double)height-(xmax-xmin);
  177.         xmin -= correct/2.;
  178.         xmax += correct/2.;
  179.         }
  180.     else
  181.         {correct=(xmax-xmin)*height/(double)width-(ymax-ymin);
  182.         ymin -= correct/2.;
  183.         ymax += correct/2.;
  184.         }
  185.     hw=height*width;
  186.     if(cycles<height) cycles=height;
  187.     if(cycles<width) cycles=width;
  188.     printf("x in (%g,%g) and y in (%g,%g)\n",xmin,xmax,ymin,ymax);
  189.     type=malloc(hw*sizeof(int));
  190.     volt=malloc((boundaries+hw)*sizeof(int)); 
  191.     if(!type || !volt)
  192.         {fprintf(stderr,"not enough workspace for %d by %d grid",height,width);
  193.         exit(1);
  194.         }
  195.     volt += boundaries;        /* allow space for the boundary potentials */
  196.     laplace(type,volt,width,height,cycles);
  197.     if(show_potentials) show_volt(type,volt,width,height);
  198. #ifdef FIXED
  199.     charges(type,volt,width,height);
  200. #endif
  201.     printf("writing result to file %s\n",outname);
  202.     if(-1==fprintf(ofile,"%10.4g%10.4g%5d    minimum and maximum x values, width\n",xmin,xmax,width)) exit(1);
  203.     if(-1==fprintf(ofile,"%10.4g%10.4g%5d    minimum and maximum y values, height\n",ymin,ymax,height)) exit(1);
  204.     yval=ymin;
  205.     for (i=0; i<height; i++)
  206.         {col=0;
  207.         for (j=0; j<width; j++)
  208.             {if(-1==fprintf(ofile,"%10.4g",volt[i*width+j]/zscale+zmin)) exit(1);
  209.             if(++col>=7)
  210.                 {if(-1==fprintf(ofile,"\n")) exit(1);
  211.                 col=0;
  212.                 }
  213.             }
  214.         if(col)
  215.             if(-1==fprintf(ofile,"\n")) exit(1);
  216.         }
  217.     fclose(ofile);
  218. }
  219.  
  220. laplace (type,volt,width,height,cycles) int *type; int *volt,width,height,cycles;
  221. {
  222. #ifdef DESMET
  223.     long start,tics();
  224. #endif
  225.     int i,j,k,hw;
  226.     hw=height*width;
  227.     if(height>width) i=height;
  228.     else i=width;
  229.     if(i>=8 && (width&1)==0 && (height&1)==0)
  230.         {laplace(type,volt,width/2,height/2,cycles/2);
  231.         i=hw; j=i/4;
  232.         while(i)
  233.             {i--; j--;
  234.             volt[i]=volt[i-1]=volt[i-width]=volt[i-width-1]=volt[j];
  235.             i--;
  236.             if(i%width==0) i -= width;
  237.             }
  238.         }
  239.     else
  240.         {for (i=0; i<hw; i++) volt[i]=4096;
  241.         }
  242.     printf("height=%d, width=%d, cycles=%d\n",height,width,cycles);
  243. #ifdef DESMET
  244.     start=tics();
  245. #endif
  246.     fill(type,volt,width,height);
  247.     k=0;
  248.     if(show_boundaries)
  249.         {for (i=height-1; i>=0; i--)
  250.             {k=i*width;
  251.             printf("%3d: ",i);
  252.             for (j=0; j<width; j++)
  253.                 printf("%3d",type[k++]);
  254.             printf("\n");
  255.             }
  256.         if(debugging) getchar();
  257.         }
  258. /*    show_volt(type,volt,width,height); */
  259.     while(cycles>-1)
  260.         {while(!csts() && cycles--)
  261.             {if(!show_resid)
  262.                 printf("\015 %d cycles to go  ",cycles);
  263.             relax(type,volt,width,height);
  264.             work += hw;
  265.             if(show_resid && ++done%5==0)
  266.                 printf("%D   %f\n",work,residual(type,volt,width,height));
  267.             }
  268.         if(cycles>-1)
  269.             {printf("\nINTERRUPTED                 ");
  270.             show_volt(type,volt,width,height);
  271.             printf("continue calculation (Y/N)? ");
  272.             if(toupper(getchar())!='Y')break;
  273.             }
  274.         }
  275. #ifdef DESMET
  276.     printf("\015 finished in %4.2f sec  \n",.01*(tics()-start));
  277. #else
  278.     printf("\015 finished               \n");
  279. #endif
  280. }
  281.  
  282. relax (type,v,w,h) int *type; int *v,w,h;
  283. {    int i,wh;
  284.     int hm1,wm1,j;
  285.     boundary_rule *rule;
  286.     long sum;
  287.  
  288. #ifdef SWITCH
  289.     i=wh=w*h;
  290.     while(i--)
  291.         {if(*type>0)
  292.             {switch(*type)
  293.                 {
  294.                 case  1: v[0]= (  v[-w] + v[-1] +     v[1] + v[w] + 2)/4; break;
  295.                 case  2: v[0]= (                      v[1] + v[w] + 1)/2; break;
  296.                 case  3: v[0]= (        2*v[-1] + 2*v[1]   + v[w] + 2)/5; break;
  297.                 case  4: v[0]= (          v[-1]            + v[w] + 1)/2; break;
  298.                 case  5: v[0]= (2*v[-w] + v[-1]          + 2*v[w] + 2)/5; break;
  299.                 case  6: v[0]= (  v[-w] + v[-1]                   + 1)/2; break;
  300.                 case  7: v[0]= (  v[-w] + 2*v[-1] + 2*v[1]        + 2)/5; break;
  301.                 case  8: v[0]= (  v[-w]         +     v[1]        + 1)/2; break;
  302.                 case  9: v[0]= (2*v[-w]         +   v[1] + 2*v[w] + 2)/5; break;
  303.                 case 10: v[0]= (2*v[ w]   + v[-1] + v[1]           + 2)/4; break;
  304.                 case 11: v[0]= (  v[-w] + 2*v[-1]          + v[ w] + 2)/4; break;
  305.                 case 12: v[0]= (            v[-1] + v[1] + 2*v[-w] + 2)/4; break;
  306.                 case 13: v[0]= (  v[-w]         + 2*v[1]   + v[ w] + 2)/4; break;
  307.                 }
  308.             }
  309.         else if (*type<0)
  310.             {rule=&bo.b_rule[-*type];
  311.             sum=0L;
  312. #ifdef DEBUG
  313. printf("\ngrid %d: ",wh-i-1); show_b(rule);
  314. show_volt(type,volt,w,h);
  315. #endif
  316.             for (j=0; j<4; j++) sum +=    volt[rule.var[j]]*(long)rule.coef[j];
  317.             v[0]=sum/SCALE;
  318. #ifdef DEBUG
  319. printf("setting potential to %5.2f\n",v[0]/zscale+zmin);
  320. #endif
  321.             }
  322.         v++; type++;
  323.         }
  324. #else
  325.     hm1=h-1;
  326.     wm1=w-1;
  327.     if(*type++)
  328.         v[0]= (                      v[1] + v[w] + 1)/2;
  329.     v++;
  330.     for (j=1; j<wm1; j++)
  331.         {if(*type++)
  332.             v[0]= (           v[-1] + v[1] + v[w] + 1)/3;
  333.         v++;
  334.         }
  335.     if(*type++)
  336.         v[0]= (           v[-1]           + v[w] + 1)/2;
  337.     v++;
  338.     for (i=1; i<hm1; i++)
  339.         {
  340.         if(*type++)
  341.             v[0]= (v[-w]            + v[1] + v[w] + 1)/3;
  342.         v++;
  343.         for (j=1; j<wm1; j++)
  344.             {if(*type++)
  345.                 v[0]= (v[-w] + v[-1] + v[1] + v[w] + 2)/4;
  346.             v++;
  347.             }
  348.         if(*type++)
  349.             v[0]= (v[-w] + v[-1]           + v[w] + 1)/3;
  350.         v++;
  351.         }
  352.     if(*type++)
  353.         v[0]= (v[-w]            + v[1]           + 1)/2;
  354.     v++;
  355.     for (j=1; j<wm1; j++)
  356.         {if(*type++)
  357.             v[0]= (v[-w] + v[-1] + v[1]           + 1)/3;
  358.         v++;
  359.         }
  360.     if(*type++)
  361.         v[0]= (v[-w] + v[-1]                     + 1)/2;
  362.     v++;
  363.  
  364. #endif
  365. }
  366.  
  367. double residual (type,volt,w,h) int *type; int *volt,w,h;
  368. {    int i,hw,n;
  369.     int hm1,wm1,j;
  370.     long r,d;
  371.  
  372.     n=0;
  373.     i=hw=w*h;
  374.     r=0L;
  375.     while(i--)
  376.         {switch(*type++)
  377.             {case 0: n++; break;
  378.             case 1: d=volt[0]-(volt[-w] + volt[-1] + volt[1] + volt[w] + 2)/4; r+=d*d; break;
  379.             case 2: d=volt[0]-(                      volt[1] + volt[w] + 1)/2; r+=d*d; break;
  380.             case 3: d=volt[0]-(           volt[-1] + volt[1] + volt[w] + 1)/3; r+=d*d; break;
  381.             case 4: d=volt[0]-(           volt[-1]           + volt[w] + 1)/2; r+=d*d; break;
  382.             case 5: d=volt[0]-(volt[-w] + volt[-1]           + volt[w] + 1)/3; r+=d*d; break;
  383.             case 6: d=volt[0]-(volt[-w] + volt[-1]                     + 1)/2; r+=d*d; break;
  384.             case 7: d=volt[0]-(volt[-w] + volt[-1] + volt[1]           + 1)/3; r+=d*d; break;
  385.             case 8: d=volt[0]-(volt[-w]            + volt[1]           + 1)/2; r+=d*d; break;
  386.             case 9: d=volt[0]-(volt[-w]            + volt[1] + volt[w] + 1)/3; r+=d*d; break;
  387.             }
  388.         volt++;
  389.         }
  390.     if (hw<=n)
  391.         return 1.;
  392.     return sqrt(((double)r)/(hw-n));
  393. }
  394.  
  395. #ifdef FIXED
  396. /*
  397.     this no longer works...could be changed to solve linear system for
  398.     the gradient at center of each grid point next to a boundary and
  399.     integrate, but it's easier to let CONTOUR calculate the integral.
  400. */
  401. charges (type,volt,width,height) int *type; int *volt,width,height;
  402. {    int v,i,plate=0,out=0,hw;
  403.     long charge;
  404.     double C,Z,q[2],plate_voltage[2],sum=0.;
  405.     hw=height*width;
  406.     printf("\nComputing surface integral of grad V over each boundary\n");
  407.     printf("(integral of -grad V dot N, where N is normal to the boundary)...\n");
  408.     printf("\nboundary    potential    surface integral     charge\n");
  409.     while(1)
  410.         {for (i=0; i<hw; i++)
  411.             if(type[i]==0)
  412.                 break;
  413.         if(i==hw) break;
  414.         v=volt[i];
  415.         plate_voltage[plate&1]=v/zscale+zmin;
  416.         charge=0.;
  417.         if(i%width) charge += v-volt[i-1];
  418.         if(i>=width) charge += v-volt[i-width];
  419.         if(i<hw-width) charge += v-volt[i+width];
  420.         if(i%width != width-1) charge += v-volt[i+1];
  421.         type[i]=23;
  422.         for (i++; i<hw; i++)
  423.             {if(type[i]==0 && v==volt[i])
  424.                 {if(i%width) charge += v-volt[i-1];
  425.                 if(i>=width) charge += v-volt[i-width];
  426.                 if(i<hw-width) charge += v-volt[i+width];
  427.                 if(i%width != width-1) charge += v-volt[i+1];
  428.                 type[i]=23;        /* don't look at this cell again */
  429.                 }
  430.             }
  431.         sum += (q[plate&1]=(charge/zscale)*PERMITTIVITY);
  432.         printf("%5d %14.2g %15.3g %15.3g\n",
  433.             plate,v/zscale+zmin,q[plate&1]/PERMITTIVITY,q[plate&1]);
  434.         plate++;
  435.         }
  436.     printf("    sum (should be zero)%12.3g %15.3g\n\n",sum/PERMITTIVITY,sum);
  437.     if(plate==2)
  438.         {q[0]=fabs(q[0]);
  439.         q[1]=fabs(q[1]);
  440.         C=MAX(q[1],q[0])/fabs(plate_voltage[1]-plate_voltage[0]);
  441.         Z=1./C/2.9979e8;
  442.         printf("characteristic impedance is %5.2f ohms\n",Z);
  443.         printf("capacitance is %7.3g farad/m\n",C);
  444.         printf("inductance is %7.3g henries/m\n",Z/2.9979e8);
  445.         }
  446. }
  447. #endif
  448. /*
  449.     the equations relating a grid point and its four neighbors, in the
  450.     absence of a boundary, is:
  451.     
  452.         ( 0  1  1  1) (Ux   )   (Uright)
  453.         ( 0 -1  1  1) (Uy   ) = (Uleft )
  454.          ( 1  0 -1  1) (Uxx  )   (Uup   )
  455.         (-1  0 -1  1) (Uhere)   (Udown )
  456.  
  457.     where    Ux is h*(d/dx U)
  458.             Uy is h*(d/dy U)
  459.             Uxx is .5*h**2*(d^2/dx^2 U))
  460. */
  461.  
  462. double lhs[16]=
  463.     {0., 1., 1., 1.,
  464.      0.,-1., 1., 1.,
  465.      1., 0.,-1., 1.,
  466.     -1., 0.,-1., 1.};
  467.  
  468. #define LEFT_SIDE(g)    ((g)%w==0)
  469. #define RIGHT_SIDE(g)    ((g)%w==w-1)
  470. #define TOP_SIDE(g)        ((g)<w)
  471. #define LOWER_SIDE(g)    ((g)>=hw-w)
  472.  
  473.     /*    fill type and voltage arrays based on user's boundary conditions */
  474. fill (type,volt,w,h) int *type; int *volt,w,h;
  475. {    boundary_rule *ru;
  476.     crossing *this,*last;
  477.     double invert(), decomp();
  478.     double a[4][4], ai[4][4], cond, condp1, *p, f;
  479.     double xt,yt,xt2,yt2;
  480.     int i,j,g,weight,zero;
  481.     int compare();
  482.     int i,k,hw,zz;
  483.  
  484.     printf("establishing boundary conditions\n");
  485.     hw=h*w;
  486.     xscale=(w-1)/(xmax-xmin);
  487.     yscale=(h-1)/(ymax-ymin);    /* same as xscale */
  488.     edge=1./xscale;
  489. /*    printf("xscale=%f, yscale=%f, edge=%f\n",xscale,yscale,edge); */
  490.     zscale=32760./5./(zmax-zmin);
  491.  
  492.     for (i=0; i<hw; i++)    type[i]=1;        /* interior */
  493.     if(grounded)
  494.         {zero=floor((0.-zmin)*zscale+.499);
  495.         for (i=w; i<hw; i+=w)    {type[i]=0; volt[i]=zero;}    /* left */
  496.         for (i=0; i<w; i++)        {type[i]=0; volt[i]=zero;}    /* top */
  497.         for (i=w-1; i<hw; i+=w)    {type[i]=0; volt[i]=zero;}    /* right */
  498.         for (i=hw-w; i<hw; i++)    {type[i]=0; volt[i]=zero;}    /* bottom */
  499.         }
  500.     else
  501.         {if(left_symmetric)
  502.             for (i=w; i<hw; i+=w)    type[i]=13;        /* left */
  503.         else
  504.             for (i=w; i<hw; i+=w)    type[i]=9;        /* left */
  505.         for (i=0; i<w; i++)        type[i]=3;            /* top */
  506.         for (i=w-1; i<hw; i+=w)    type[i]=5;            /* right */
  507.         if(lower_symmetric)
  508.             for (i=hw-w; i<hw; i++)    type[i]=12;        /* bottom */
  509.         else
  510.             for (i=hw-w; i<hw; i++)    type[i]=7;        /* bottom */
  511.         type[0]=2; type[w-1]=4; type[hw-w]=8; type[hw-1]=6;        /* corners */
  512.         }
  513.     i=k=0;
  514.     lc=2;    /* bo.cross_list is empty */
  515.     while(i<n)
  516.         {xt=x[i]-xmin;
  517.         yt=y[i]-ymin;
  518.         zz=floor((p_voltage[k]-zmin)*zscale+.499);
  519.         volt[-k-1]=zz;
  520.         do    {xt2=x[i]-xmin;
  521.             yt2=y[i]-ymin;
  522.             mark(w,h,xt,yt,xt2,yt2,-k-1);
  523. /*
  524.             printf("%f,%f -> ",(x[i]-xmin)*xscale,(y[i]-ymin)*yscale);
  525.             printf("mark(%d,%d, %d,%d, %d,%d, %d)\n",w,h,xx,yy,xx2,yy2,zz); 
  526. */
  527.             xt=xt2; yt=yt2; 
  528.             i++;
  529.             } while(i<=p_data[k]);
  530.         k++;
  531.         }
  532. /*    printf("the first of the %d crossings before sorting\n",lc-1);
  533.     for (i=2; i<9; i++) show_c(bo.cross_list[i]); getchar();    */
  534. /*----------------  sort list  -------------*/
  535.     hsort(&bo.cross_list[2],lc-2,sizeof(crossing),&compare);
  536. /*----------------  solve each problem  -------------*/
  537. #ifdef DEBUG
  538.     printf("crossings after sorting\n");
  539.     for (i=2; i<9; i++) show_c(bo.cross_list[i]); getchar();
  540. #endif
  541.     lb=1;    /* bo.b_rule is empty */
  542.     this=&bo.cross_list[2];
  543.     last=&bo.cross_list[lc];
  544.     while(this<last)
  545.         {p=a;
  546.         for (i=0; i<16; i++) p[i]=lhs[i];
  547.         g=this->grid;
  548. #ifdef DEBUG
  549. printf("rule %d, boundary conditions for grid point %d\n",lb,g);
  550. #endif
  551.         volt[g]=floor((p_voltage[-this->volt-1]-zmin)*zscale+.499);
  552.         ru=&bo.b_rule[lb];
  553.         ru->var[0]=g+1;
  554.         ru->var[1]=g-1;
  555.         ru->var[2]=g+w;
  556.         ru->var[3]=g-w;
  557.         do    {j=this->direct;
  558. #ifdef DEBUG
  559. printf("        "); show_c(this);
  560. #endif
  561.             f=this->dist;
  562.             ru->var[j]=this->volt;
  563.             a[j][0] *= f;
  564.             a[j][1] *= f;
  565.             a[j][2] *= f*f;
  566.             this++;
  567.             } while (this->grid == g);
  568. /*        show_m("boundary condition matrix",a); getchar();    */
  569.         cond=invert(4,4,a,ai);
  570. /*        printf("condition number is %f\n",cond);
  571.         show_m("matrix inverse",ai);    */
  572.         condp1=cond+1.;
  573.         if(cond==condp1)
  574.             {fprintf(stderr,
  575.                 "can\'t solve linear system for boundary conditions\n");
  576.             exit(1);
  577.             }
  578.         for (i=0; i<4; i++) ru->coef[i]=ai[3][i]*SCALE+.5;
  579.         if(LOWER_SIDE(g))
  580.             {if(lower_symmetric)
  581.                 ru->coef[2] += ru->coef[3];
  582.             else
  583.                 {weight=ru->coef[2]=(ru->coef[2] + ru->coef[3])/4;
  584.                 weight = SCALE - weight + ru->coef[0]+ru->coef[1];
  585.                 ru->coef[0] += weight/2;
  586.                 ru->coef[1] += weight - weight/2;
  587.                 }
  588.             ru->coef[3]=0;
  589.             }
  590.         if(LEFT_SIDE(g))
  591.             {if(left_symmetric)
  592.                 ru->coef[0] += ru->coef[1];
  593.             else
  594.                 {weight=ru->coef[0]=(ru->coef[0] + ru->coef[1])/4;
  595.                 weight = SCALE - weight + ru->coef[2] + ru->coef[3];
  596.                 ru->coef[2] += weight/2;
  597.                 ru->coef[3] += weight - weight/2;
  598.                 }
  599.             ru->coef[1]=0;
  600.             }
  601.         if(TOP_SIDE(g))
  602.             {weight=ru->coef[3]=(ru->coef[2] + ru->coef[3])/4;
  603.             weight = SCALE - weight + ru->coef[0]+ru->coef[1];
  604.             ru->coef[1] += weight/2;
  605.             ru->coef[0] += weight - weight/2;
  606.             ru->coef[2]=0;
  607.             }
  608.         if(RIGHT_SIDE(g))
  609.             {weight=ru->coef[1]=(ru->coef[1] + ru->coef[0])/4;
  610.             weight = SCALE - weight + ru->coef[2]+ru->coef[3];
  611.             ru->coef[2] += weight/2;
  612.             ru->coef[3] += weight - weight/2;
  613.             ru->coef[0]=0;
  614.             }
  615. #ifdef DEBUG
  616. printf("        result is ");show_b(ru);
  617. #endif
  618.         type[g]=-lb;
  619.         lb++;
  620.         }
  621.  
  622. }
  623.  
  624. compare(a,b) crossing *a,*b;
  625. {    return (a->grid - b->grid);
  626. }
  627.  
  628. mark (width,height,x1,y1,x2,y2,v) int width,height,v; double x1,y1,x2,y2;
  629. {    int q,dx,dy,ax,ay,decide,up,over,i;
  630.     int nx,nx2,ny,ny2;
  631.     double del,xe,ye,f,t;
  632. /*
  633.     if(x1<0||x1>=width||x2<0||x2>=width||y1<0||y1>=height||y2<0||y2>=height)
  634.         {fprintf(stderr,"calling sequence error: mark(%d, %d,%d, %d,%d, %d)\n",
  635.             width,x1,y1,x2,y2,v);
  636.         fprintf(stderr,"(point outside 0<=x<%d or 0<=y<%d)\n",width,height);
  637.         exit(1);
  638.         }
  639. */
  640. #ifdef DEBUG
  641. printf("line from  %5.2f,%5.2f  to  %5.2f,%5.2f at potential %5.2f\n",
  642. x1,y1,x2,y2,p_voltage[-v-1]);
  643. #endif
  644.     if(x2<x1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
  645.     nx2=x2*xscale;
  646.     del=x2-x1;
  647.     for (nx=x1*xscale+1; nx<=nx2; nx++)
  648.         {xe=nx*edge;
  649.         ye=(y1*(x2-xe) + y2*(xe-x1))/del;
  650.         ny=ye*yscale;
  651.         f=ye/edge-ny;
  652. #ifdef DEBUG
  653. printf("     crosses y grid line at  %f,%f  ->  %d,%d at distance %f\n",xe,ye,nx,ny,f);
  654. #endif
  655.         register_crossing(width,nx,ny,f,2,v);
  656.         register_crossing(width,nx,ny+1,1.-f,3,v);
  657.         }
  658.     if(y2<y1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
  659.     ny2=y2*yscale;
  660.     del=y2-y1;
  661.     for (ny=y1*yscale+1; ny<=ny2; ny++)
  662.         {ye=ny*edge;
  663.         xe=(x1*(y2-ye) + x2*(ye-y1))/del;
  664.         nx=xe*xscale;
  665.         f=xe/edge-nx;
  666. #ifdef DEBUG
  667. printf("     crosses x grid line at  %f,%f  ->  %d,%d at distance %f\n",xe,ye,nx,ny,f);
  668. #endif
  669.         register_crossing(width,nx,ny,f,0,v);
  670.         register_crossing(width,nx+1,ny,1.-f,1,v);
  671.         }
  672. }
  673.  
  674. register_crossing(width,nx,ny,f,direction,v)
  675. int width,nx,ny,direction,v;
  676. double f;
  677. {    crossing *new;
  678.  
  679.     if(lc>=MAX_CROSS) {fprintf(stderr,"too many boundary crosspoints\n"); exit(1);}
  680.     new=&bo.cross_list[lc++];
  681.     new->grid=nx+ny*width;
  682.     new->direct=direction;
  683.     new->dist=f;
  684.     new->volt=v;
  685. }
  686.  
  687. show_volt (type,volt,w,h) int *type; int *volt,w,h;
  688. {    int i,j,k;
  689.     k=0;
  690. printf("boundaries...");
  691. for (i=1; i<=boundaries; i++) printf("%d: %5.2f   ",i,volt[-i]/zscale+zmin);
  692. printf("\n");
  693. /*    printf("(press any key to abort printout)\n");
  694.     for (i=0; i<h && !csts(); i++)
  695. */
  696.     for (i=height-1; i>=0; i--)
  697.         {k=i*w;
  698.         printf("%3d: ",i);
  699.         for (j=0; j<w; j++)
  700.             {printf("%6.2f",volt[k++]/zscale+zmin);
  701.             /* printf("%6d",volt[k++]); */
  702.             }
  703.         printf("\n");
  704.         }
  705.     for (i=height-1; i>=0; i--)
  706.         {k=i*w+j;
  707.         printf("%3d: ",i);
  708.         for (j=0; j<width; j++)
  709.             printf("%3d",type[k++]);
  710.         printf("\n");
  711.         }
  712.     if(debugging) getchar();
  713. }
  714.  
  715. read_data(argc,argv) int argc; char **argv;
  716. {    int i,j,length;
  717.     double xx,yy,d,*pd,sum;
  718.     char *s,*t;
  719. #define BUFSIZE 200
  720.     static char buf[BUFSIZE];
  721.  
  722.     x=malloc(ENTRIES*sizeof(double));
  723.     y=malloc(ENTRIES*sizeof(double));
  724.     if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
  725.     if(argc>1 && streq(argv[1],"?")) help();
  726.     if(argc<=1 || *argv[1]=='-') file=stdin;
  727.     else
  728.         {if(argc>1)
  729.             {file=fopen(argv[1],"r");
  730.             if(file==0) {printf("file %s not found\n",argv[1]); exit();}
  731.             argc--; argv++;
  732.             }
  733.         else help();
  734.         }
  735.     argc--; argv++;
  736.     while(argc>0)
  737.         {i=get_parameter(argc,argv);
  738.         argv=argv+i; argc=argc-i;
  739.         }
  740.     p_data[0]=-1;
  741.     if(height<2||width<2||cycles<2)
  742.         {fprintf(stderr,"parameter too small: height=%d, width=%d, cycles=%d",
  743.             height,width,cycles);
  744.         exit(1);
  745.         }
  746.     i=0;
  747.     while(i<ENTRIES)
  748.         {if(fgets(buf,BUFSIZE,file)==0) break;
  749.         t=buf;
  750.         while(*t && isspace(*t)) t++;
  751.         if(*t == 0) continue;        /* skip blank lines */
  752.         buf[strlen(buf)-1]=0;        /* zero out the line feed */
  753.         if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
  754.         sscanf(buf,"%F %F",&x[i],&y[i]);
  755.         if (x[i]<xmin) xmin=x[i];
  756.         if (x[i]>xmax) xmax=x[i];
  757.         if (y[i]<ymin) ymin=y[i];
  758.         if (y[i]>ymax) ymax=y[i];
  759.         s=buf;                      /* start looking for label */
  760.         while(*s==' ')s++;            /* skip first number */
  761.         while(*s && (*s!=' '))s++;
  762.         while(*s==' ')s++;            /* skip second number */
  763.         while(*s && (*s!=' '))s++;
  764.         while(*s==' ')s++;
  765.         if((length=strlen(s))&&(boundaries<MAXBOUNDARIES))
  766.             {if(*s=='\"') s++;           /* label in quotes */
  767.             sscanf(s,"%F",&p_voltage[boundaries]);
  768.             if(p_voltage[boundaries]<zmin) zmin=p_voltage[boundaries];
  769.             if(p_voltage[boundaries]>zmax) zmax=p_voltage[boundaries];
  770.             p_data[boundaries++]=i;
  771.             }
  772.         i++;
  773.         }
  774.     if(grounded)
  775.         {if(0.<zmin) zmin=0.;
  776.         if(0.>zmax) zmax=0.;
  777.         }
  778.     n=i;
  779.     if(!boundaries || p_data[boundaries-1]!=n-1)
  780.         {p_data[boundaries]=i-1;
  781.         p_voltage[boundaries++]=0.;
  782.         }
  783.     xmin -= margin;
  784.     xmax += margin;
  785.     ymin -= margin;
  786.     ymax += margin;
  787. }
  788.  
  789.  
  790. /* get_parameter - process one command line option
  791.         (return # parameters used) */
  792. get_parameter(argc,argv) int argc; char **argv;
  793. {    int i;
  794.     double temp;
  795.  
  796.     if(streq(*argv,"-d")) {debugging=1; return 1;}
  797.     else if(streq(*argv,"-h"))
  798.         {if((argc>1) && numeric(argv[1])) height=atoi(argv[1]);
  799.         return 2;
  800.         }
  801.     else if(streq(*argv,"-w"))
  802.         {if((argc>1) && numeric(argv[1])) width=atoi(argv[1]);
  803.         return 2;
  804.         }
  805.     else if(streq(*argv,"-g"))
  806.         {grounded++;
  807.         return 1;
  808.         }
  809.     else if(streq(*argv,"-m"))
  810.         {i=get_double(argc,argv,1,&margin,&margin,&margin);
  811.         if(margin<0.) margin=0.;
  812.         return i;
  813.         }
  814.     else if(streq(*argv,"-n"))
  815.         {if((argc>1) && numeric(argv[1])) cycles=atoi(argv[1]);
  816.         return 2;
  817.         }
  818.     else if(streq(*argv,"-r"))
  819.         {show_resid++;
  820.         return 1;
  821.         }
  822.     else if(streq(*argv,"-b"))
  823.         {show_boundaries++;
  824.         return 1;
  825.         }
  826.     else if(streq(*argv,"-p"))
  827.         {show_potentials++;
  828.         return 1;
  829.         }
  830. /*    else if(streq(*argv,"-a"))
  831.         {i=get_double(argc,argv,1,&temp,&temp,&temp);
  832.         alpha=temp;
  833.         alphap1=alpha+1;
  834.         return i;
  835.         }
  836. */
  837.     else if(streq(*argv,"-x"))
  838.         {i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
  839.         x_arguments=i-1;
  840.         return i;
  841.         }
  842.     else if(streq(*argv,"-y"))
  843.         {i=get_double(argc,argv,2,&ymin,&ymax,&ymax);
  844.         y_arguments=i-1;
  845.         return i;
  846.         }
  847.     else if(streq(*argv,"-xs"))
  848.         {left_symmetric++;
  849.         return 1;
  850.         }
  851.     else if(streq(*argv,"-ys"))
  852.         {lower_symmetric++;
  853.         return 1;
  854.         }
  855.     else gripe(argv);
  856. }
  857.  
  858. get_double(argc,argv,permitted,a,b,c)
  859. int argc,permitted; char **argv; double *a,*b,*c;
  860. {    int i=1;
  861.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
  862.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
  863.     if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
  864.     return i;
  865. }
  866.  
  867. int streq(a,b) char *a,*b;
  868. {    while(*a)
  869.         {if(*a!=*b)return 0;
  870.         a++; b++;
  871.         }
  872.     return 1;
  873. }
  874.  
  875. gripe_arg(s) char *s;
  876. {    fprintf(stderr,"argument missing for switch %s",s);
  877.     help();
  878. }
  879.  
  880. gripe(argv) char **argv;
  881. {    fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
  882.     help();
  883. }
  884.  
  885. help()
  886. {    fprintf(stderr,"laplace   version %s",VERSION);
  887.     fprintf(stderr,"\nusage: laplace [file][options]\n");
  888.     fprintf(stderr,"options are:\n");
  889. /*    fprintf(stderr,"  -a  num            acceleration factor (default %d) \n",alpha); */
  890.     fprintf(stderr,"  -b               display boundaries\n");
  891.     fprintf(stderr,"  -g               edge of grid is grounded (potential 0)\n");
  892.     fprintf(stderr,"  -p               display calculated potentials\n");
  893.     fprintf(stderr,"  -r               calculate and display residuals\n");
  894.     fprintf(stderr,"  -m  margin       amount to expand grid beyond given boundaries (default 0) \n");
  895.     fprintf(stderr,"  -n  num          number of relaxation cycles (default %d) \n",cycles);
  896.     fprintf(stderr,"  -h  num          height of grid in cells (default %d)\n",height);
  897.     fprintf(stderr,"  -w  num          width of grid in cells (default %d)\n",width);
  898.     fprintf(stderr,"  -x  [min [max]]  specify x range \n");
  899.     fprintf(stderr,"  -y  [min [max]]  specify y range \n");
  900.     fprintf(stderr,"  -xs              symmetric across the line x=xmin \n");
  901.     fprintf(stderr,"  -ys              symmetric across the line y=ymin \n");
  902.     exit();
  903. }
  904.  
  905. numeric(s) char *s;
  906. {    char c;
  907.     while(c=*s++)
  908.         {if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
  909.         return 0;
  910.         }
  911.     return 1;
  912. }
  913.  
  914. #ifdef DESMET
  915.  
  916. /*
  917.     tics - report hundreths of a second since midnight
  918. */
  919.  
  920. long tics()
  921. {    int hr,min,sec,hun;
  922.  
  923.     _timer(&hr,&min,&sec,&hun);
  924.     return (( (long) hr*60+min)*60+sec)*100+hun;
  925. }
  926.  
  927. _timer()
  928. {
  929. #asm
  930.     mov    ah,2ch
  931.     int    21h
  932.     mov    bx,[bp+4]
  933.     mov    [bx],ch        ;hours
  934.     mov    byte [bx+1],0
  935.     mov    bx,[bp+6]
  936.     mov    [bx],cl        ;minutes
  937.     mov    byte [bx+1],0
  938.     mov    bx,[bp+8]
  939.     mov    [bx],dh        ;seconds
  940.     mov    byte [bx+1],0
  941.     mov    bx,[bp+10]
  942.     mov    [bx],dl        ;hundreths
  943.     mov    byte [bx+1],0
  944. #endasm
  945. }
  946.  
  947. #endif
  948.  
  949. show_m(s,a) char *s; double *a;
  950. {    int i,j;
  951.     printf("%s\n",s);
  952.     for (i=0; i<4; i++)
  953.         {for (j=0; j<4; j++) printf("%8.2f ",*a++);
  954.         printf("\n");
  955.         }
  956. }
  957.